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ABSTRACT 


A comparison  of  numerical  methods  utilized  by  the  finite 
element  technique  for  solving  a nonlinear  nuclear  reactor 
dynamics  problem  was  conducted.  Using  the  Crank-Nicolson, 
DVOGER  (Gear)  and  Implicit  Gear  methods,  the  results  showed 
the  Implicit  to  be  the  superior  method  investigated.  This 
is  based  on  the  fact  that  all  three  methods  yielded  the  same 
steady  state  solutions;  but,  the  Implicit  Gear  method  used 
significantly  less  CPU  time  and  comparable  storage  to  Crank- 
Nicolson.  This  was  particularly  apparent  as  the  degrees  of 
freedom  were  increased.  In  addition,  the  transient  solution 
in  all  cases  was  better  than  that  obtained  in  Crank-Nicolson 
and  compared  favorably  to  that  of  Gear's  method. 

The  other  noteworthy  result  was  in  the  effect  of  the 

error  criterion  on  solution.  It  was  shown  that  for  a range 

- 4 

of  error  from  10  to  1.0,  the  steady  state  solution  value 
remained  the  same.  This  results  in  a significant  reduction 
in  computer  processing  time  since  the  time  required  decreases 
substantially  as  the  error  conditions  imposed  are  relaxed. 
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I.  INTRODUCTION 


A.  PURPOSE 

This  research  project  has  been  undertaken  to  compare 
several  numerical  methods  of  solving  a nonlinear  nuclear 
reactor  dynamics  problem.  Three  methods  have  been  investi- 
gated in  this  thesis,  including  the  Crank-Nicolson , the 
DVOGER  (Gear)  and  the  Implicit  Gear  methods  of  solution. 

A nuclear  reactor  dynamics  problem  with  temperature- 
dependent  feedback,  when  it  entails  either  a non-homogeneous 
or  multi-region  reactor,  results  in  a nonlinear  field  equa- 
tion in  space  and  time.  This  problem  does  not  lend  itself 
to  solution  by  analytical  means  [5,  6],  However,  when  the 
physical  and  neutronic  properties  of  the  problem  are  known, 
a model  can  be  formulated  using  the  finite  element  method 
which  will  yield  the  transient  and  steady  state  flux  solu- 
tions. In  particular,  the  partial  differential  equations 
investigated  were  of  the  form 

aTt  = bv2  * + CIP  -v* 2 Cl) 

where  a,  b,  c,  and  w are  constants  and  ip  (r,  z,  t)  is  the 
flux.  The  finite  element  method  reduces  Equation  (1)  to 
the  system  of  ordinary  differential  equations 
N . N 

j=l  Aij  Yj(t)  “ j£i  Bij  i|>jU)  i - 1,  ...»  N (2) 
when  the  nonlinear  term  is  linearized.  Nguyen  and  Salinas  [5] 
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and  more  recently  Olsen  [6]  discuss  the  problem  and  methods 
of  solution. 

In  this  work,  a comparison  of  three  numerical  integration 
schemes  for  Equation  (2)  was  made.  The  comparisons  will  be 
made  on  several  bases:  computer  storage  requirements,  com- 

puter processing  time,  rapidity  with  which  solution  was 
obtained  and  the  relative  accuracy  of  the  solution. 

At  steady  state,  it  is  expected  that  all  three  methods 
will  provide  the  same  solution  value.  This  is  due  to  the 
fact  that  at  steady  state  the  time  derivative  of  the  flux 

t 

( ) , Equation  (2)  , is  zero. 

In  order  to  explore  the  relative  value  of  the  various 
methods,  the  model  has  been  discretized  into  finite  element 
grids  of  various  degrees  of  freedom  (DOF) . The  effect  on 
the  solution  by  finer  discretizations  of  the  finite  element 
has  been  investigated  to  determine  if  the  various  methods  of 
solution  are  similarly  influenced  to  provide  a better  solution 
for  a finer  mesh.  To  have  a method  of  solution  relatively 
independent  of  mesh  size  would  greatly  reduce  computer  storage 
and  processing  time  requirements  since  a larger  grid  with 
fewer  elements  could  be  utilized. 

In  order  to  test  the  flexibility  of  the  various  equation 
solvers,  an  initial  disturbance  was  introduced  at  different 
points  in  the  model.  This  provided  both  a check  of  the 
ability  of  the  method  to  accept  a random  disturbance  as  well 
as  information  on  how  rapidly  a solution  is  obtained  with  a 
particular  disturbance  input.  Two  nodal  points  of  the  system 


were  considered,  a point  at  the  origin  and  a point  on  the 
core-reflector  interface.  This  was  done  to  determine  if  the 
tracking  ability  was  consistent  throughout  the  model. 

Additional  comparisons  investigated  include  the  effects 
of  the  convergence  criterion  on  the  solution  for  all  methods 
and  the  effects  of  the  size  of  the  time  step  on  the  Crank- 
Nicolson  method. 

Modifications  have  been  made  to  the  programs  provided  in 
Ref.  [6]  which  was  the  basis  of  this  research  project.  In 
that  work,  core  properties  were  arbitrarily  assigned  to  the 
reflector  elements  at  the  interface.  This  occurred  because 
the  material  and  nuclear  properties  were  provided  at  the  nodal 
points  rather  than  by  elements.  The  properties  were  intro- 
duced in  this  manner  as  a means  of  computer  storage  reduction. 
Regardless  of  the  mesh  size,  there  were  always  fewer  nodes 
than  elements.  In  this  project,  all  properties  were  provided 
on  an  element  basis  to  eliminate  this  discrepancy  while  at 

; 

the  same  time  sacrificing  the  minor  storage  savings  it  repre- 
sented. 


II.  DATA  GENERATORS 


A.  GENERAL 

The  data  generators  formulated  in  this  work  are  useable 
only  for  regular  rectangular  grids  (Figures  1-3).  Having 
selected  the  number  of  horizontal  and  vertical  points  for 
the  discretization  of  the  model,  the  data  generators  will 
provide  the  numbering  of  each  nodal  point,  the  horizontal 
and  vertical  position  of  each  node  and  the  nodal  neighbors 
of  each  node. 

In  addition,  the  outer  boundary  nodes,  at  which  the 
dynamic  flux  is  zero,  will  be  numbered  last  in  order  to 
reduce  the  number  of  equations  required  by  the  particular 
method  of  solution  being  used.  This  wxll  substantially 
decrease  the  storage  requirements.  For  example,  in  a one 
hundred  thirty- two  node  discretization,  only  one  hundred  ten 
equations  will  be  solved. 

B.  PROPERTY  INPUTS 

1.  Purpose 

A simple  data  generator  has  been  provided  which  will 
produce  a data  deck  containing  the  physical  properties  of 
the  reactor  core  and  reflector  in  the  format  required  by  the 
Crank-Nicolson  and  DVOGER  (Gear)  methods.  These  properties 
are  neutron  velocity,  neutron  diffusion  length,  neutron 
absorption  cross-section,  neutron  fission  cross-section  and 
reactivity  temperature  coefficient. 
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This  generator  will  accommodate  any  size  mesh  and 
will  process  an  indefinite  number  of  grids  simultaneously. 

The  user  must  provide  this  routine  with  a data  deck 
which  contains  the  total  number  of  elements  in  the  grid  and 
the  type  of  each  element,  core  or  reflector.  This  has  been 
simplified  by  the  designation  of  all  core  elements  by  ITYPE=0 
and  reflector  elements  by  ITYPE=1. 

2 . Programming  Notes 

a.  The  property  values  in  the  two  regions  must  be 
provided  as  an  integral  part  of  the  program. 

b.  The  type  of  element,  core  or  reflector,  must  be 
provided. 

c.  The  routine  will  process  an  unlimited  number  of 
models.  It,  therefore,  must  be  halted  when  all  data  has  been 
utilized.  This  is  done  by  specifying  the  final  value  of  the 
number  of  elements  in  an  IF-STOP  statement. 

3.  Parameters 

a.  NEL  - number  of  elements  in  the  discretized  model 

b.  ITYPE  - type  of  element;  ITYPE*0  is  a core  element 
and  ITYPE=1  is  a reflector  element. 

c.  D - vector  containing  the  diffusion  length  of  the 
core  and  reflector  elements. 

d.  SGA  - vector  containing  the  absorption  cross- 
section  of  the  core  and  reflector  elements. 

e.  SGF  - vector  containing  the  fission  cross-section 
of  the  core  and  reflector  elements. 

f.  V - vector  containing  the  neutron  velocity. 
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g.  ALPHA  - vector  containing  the  reactivity 
temperature  coefficient. 


C.  NODAL  POINT  COORDINATES  AND  ELEMENT 
NODAL  POINT  CONNECTIVITY 

1 . Purpose 

Nodal  point  positioning  in  the  model  is  readily 
obtained  in  data  deck  form  from  this  routine.  Once  the  model 
to  be  investigated  has  been  discretized  into  the  finite  ele- 
ment grid  desired,  the  user  provides  the  number  of  vertical 
and  horizontal  nodal  points  and  their  dimensional  position 
along  the  axes.  By  the  use  of  a nested  loop,  the  nodal  points 
are  consecutively  numbered  and  assigned  the  proper  coordinate 
dimensions . 

The  element  connectivity,  that  is,  the  nodal  points 
for  each  element,  is  also  resolved  by  the  use  of  a nested 
loop  and  several  counters.  Each  iteration  will  yield  two 
elements  with  their  respective  nodal  point  boundaries.  This 
will  continue  until  the  nodal  points  for  each  element  have 
been  computed. 

The  output  of  this  routine  will  consist  of  two  data 
decks.  The  first  will  provide  the  nodal  point  with  its  ver- 
tical and  horizontal  coordinates.  The  second  will  yield  the 
element  nodal  point  connectivity. 

Omitted  from  the  element  connectivity  data  deck,  but 
required  by  the  thesis  program,  is  the  type  of  element,  core 
or  reflector.  This  must  be  added  to  the  deck  individually. 
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Example  of  output  (See  Figure  1) : 


nodal  point 
1 


nodal  point  coordinates 
horizontal  position  vertical  position 
0.0  0.0 


element 

1 


element  connectivity 
nodal  point  connectivity 
19  10 


type  element* 
0 


*Note:  must  be  added  manually  to  data  card. 

2.  Programming  Notes 

a.  The  maximum  number  of  vertical  and  horizontal 
nodal  points  in  the  discretized  model  must  be  provided. 

b.  The  horizontal  and  vertical  positions  of  the 
discretization  must  be  provided. 

c.  The  boundary  points  of  the  model  will  be  numbered 
such  that  these  nodal  points  are  the  last  in  the  numerical 
sequence . 

d.  The  routine  is  written  such  that  it  will  process 
an  indefinite  number  of  models.  Therefore,  it  must  be  halted 
when  all  the  data  has  been  utilized.  This  will  be  accomplished 
by  specifying  the  maximum  number  of  vertical  points  in  the 
model  in  an  IF-STOP  statement. 

3.  Parameters 

a.  NH  - maximum  number  of  horizontal  nodal  points 
in  the  discretized  model. 

b.  NV  - maximum  number  of  vertical  nodal  points  in 
the  discretized  model. 
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c.  X - vector  containing  the  horizontal  positions 
of  the  NH  points. 

d.  Y - vector  containing  the  vertical  positions  of 
the  NV  points. 

e.  SYSNOD  - vector  providing  the  total  number  of 
nodal  points. 

f.  R - vector  providing  the  radial  position  of  the 
nodal  points. 

g.  Z - vector  providing  the  vertical  position  of 
the  nodal  points. 

h.  ELNOD  - array  providing  the  nodal  point  boundaries 
for  each  element. 

i.  NEL  - the  total  number  of  elements  in  the  discre- 
tized model. 

D.  NODAL  NEIGHBOR  CONNECTIVITY 
1.  Purpose 

This  routine  produces  a data  deck  for  the  Crank- 
Nicolson  and  Implicit  Gear  methods  containing  each  nodal 
point  and  the  nodal  points  connected  to  it  by  the  finite 
element  discretization. 

The  regular  rectangular  grids  are  such  that  no  more 
than  six  nodal  points  contribute  to  any  one;  therefore,  an 
Nx7  array  can  be  formulated  for  nodal  neighbor  connectivity. 
For  example,  nodal  point  26  in  the  112  element  grid  would  be 
stored  as  follows  (Figure  4) : 

26  17  18  27  35  34  25 
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and  nodal  point  58 


58  49  50  59  57  0 0. 

As  written,  zeros  are  inputted  whenever  a nodal  point  does 
not  have  six  neighbors  to  maintain  the  symmetry  of  the  matrix. 
This  could,  however,  be  compacted  further  by  reading  the  data 
as  a single  vector  with  a starting  point  counter  to  indicate 
when  the  next  nodal  point  has  been  reached  in  the  sequence. 
This  would  eliminate  the  need  for  the  zeros  to  be  supplied. 

This  routine  will  provide  a data  deck  for  any 
rectangular  grid  and  for  an  indefinite  number  of  grids. 

2 . Algorithm 

This  algorithm  provides  a dense,  compacted  matrix 
which  reduces  the  storage  size  required  by  the  main  program. 
The  use  of  a finite  element  method  allows  a certain  amount  of 
compacting,  i.e.,  the  banded  matrix.  This  is  a function  of 
the  finite  element  size,  that  is,  the  greater  the  number  of 
nodal  points,  the  larger  the  band.  The  size  of  the  band  is 
determined  by  the  largest  difference  between  nodal  neighbors 
plus  one.  For  example,  in  Figure  1,  node  11  has  a maximum 
difference  of  20  - 2 ■ 18,  while  at  node  16  the  difference 
is  43  - 7 = 36  which  is  the  largest  difference  in  the  forty- 
five  node  system.  As  the  number  of  nodal  points  increases, 
this  maximum  band  width  also  increases.  In  Figure  2,  the 
maximum  difference  is  at  node  16  (67  - 7 = 60) ; and,  in 
Figure  3,  the  maximum  difference  of  124  - 10  * 114  exists 
at  node  22.  Therefore,  although  the  size  of  the  system  is 
reduced  from  NxN  to  Nxq,  the  size  is  not  constant  and  may 


19 


not  be  small.  In  general,  if  one  used  banded  storage,  the 
numbering  would  proceed  consecutively  in  the  direction  of 
fewest  nodes.  However,  this  would  require  the  identification 
of  those  nodes  on  the  boundary  since  they  are  of  constant 
value  and  do  not  require  integration. 

In  the  particular  scheme  of  discretization  utilized 
in  this  project,  no  nodal  point  has  more  than  six  nodal  point 
contributors.  This  allows  the  matrix  to  be  compacted  further 
from  Nxq  to  Nx7. 

3 . Programming  Notes 

a.  The  total  number  of  nodal  points,  the  maximum 
number  of  nodal  point  contributors  and  the  maximum  number  of 
horizontal  and  vertical  nodal  points  must  be  provided. 

b.  The  routine  will  satisfy  any  rectangular  grid. 
However,  it  must  be  instructed  when  a sufficient  number  of 
nodal  points  have  been  generated.  This  is  accomplished  by 
using  the  total  number  of  nodal  points  desired  in  an  IF-GO  TO 
statement. 

c.  This  routine  requires  a positive  means  of  stopping. 
This  is  accomplished  by  the  use  of  the  maximum  number  of 
vertical  nodal  points  in  the  last  data  set  in  an  IF-STOP 
statement . 

d.  The  output  of  this  routine  will  place  the  central 
nodal  point  first  with  the  contributors  following. 

4 . Parameters 

a.  NV  - number  of  vertical  nodal  points  in  the  model. 

b.  NH  - number  of  horizontal  nodal  points  in  the  model 


c.  NVH  - total  number  of  nodal  points  in  the  model. 


NVH  3 NV  x NH. 

d.  LCON  - maximum  number  of  contributing  nodal  points. 
(LCON  - 7) . 

e.  MNOD  - an  array,  NVH  x LCON,  providing  the  central 
and  contributing  nodal  points. 
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III.  CRANK-NICOLSON  METHOD  OF  SOLUTION 


A.  DESCRIPTION 

The  Crank-Nicolson  formulation  is  a numerical  method 
originally  presented  by  J.  Crank  and  P.  Nicolson  in  1947. 
Crank-Nicolson,  being  an  implicit  method,  does  not  require 
the  inverse  of  the  matrix  to  be  calculated.  Therefore, 
advantage  is  taken  of  the  sparse  matrix  inherently  provided 
by  the  finite  element  method.  If  the  inverse  of  the  sparse 
matrix  is  formed,  it  will  be  full. 

The  Crank-Nicolson  method  will  solve  the  set  of  linear 
differential  equations  represented  by 
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Crank-Nicolson,  after  some  algebra,  yields 
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The  implicit  system  (4)  may  be  solved  by  an  iterative  process; 
in  this  case,  the  Gauss-Seidel  method  is  used  to  solve  the 
set  of  simultaneous  algebraic  equations  (4)  formulated.  Of 
all  the  stable  numerical  methods  in  the  case  of  single  step 
implicit,  Crank-Nicolson  has  the  smallest  truncation  error  [2] . 
The  Gauss-Seidel  method  of  iteration  is  continuously  using  the 
newest  solution  values  allowing  convergence  to  the  solution 
to  be  the  most  rapid. 
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B.  EFFECT  OF  DEGREES  OF  FREEDOM 

As  might  be  expected,  as  the  grid  size  became  finer,  the 
computer  core  requirements  increase.  In  addition,  more  time 
is  required  for  the  problem  to  reach  a steady  state  solution. 
At  the  same  time,  as  shown  in  Figures  5-7,  the  solution  values 
converged,  as  the  mesh  size  became  finer,  to  a more  accurate 
solution. 

As  shown  in  Table  I and  Figure  8,  storage  and  processing 
time,  that  is,  CPU  time  to  reach  a steady  state  solution,  in- 
creased markedly  between  a forty- five  node  grid  and  a one 
hundred  thirty-two  node  grid.  This  yielded  a fifty-nine 
percent  increase  in  storage  and  better  than  a five  hundred 
percent  increase  in  processing  time.  At  the  same  time  there 
is  an  improvement  in  the  solution  of  only  four  and  one-half 
percent . 

Comparison  to  a seventy-two  node  grid  yielded  far  more 
satisfying  results.  At  seventy- two  nodes,  there  is  an  eighty 
percent  increase  in  CPU  time  and  a sixteen  percent  increase 
in  storage  requirements  while  obtaining  a two  and  one-half 
percent  improvement  in  the  solution. 

Similar  results  were  obtained  in  problem  time  to  solution. 
For  one  hundred  thirty-two  nodes,  the  time  more  than  doubles; 
while  for  seventy- two  nodes,  the  time  was  increased  by  twenty 
percent . 

C.  EFFECT  OF  ERROR  CRITERION 

The  error  criterion  used  throughout  this  thesis  informs 
the  integration  routine  when  it  has  achieved  sufficient 
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agreement  between  iterates  (i.e.,  when  it  can  stop  and  go 
to  the  next  iteration  step) . 

In  Crank-Nicolson,  the  error  criterion  simply  consisted 
of  the  difference  in  successive  Gauss-Seidel  iterates  divided 
by  the  current  iterative  value. 

- 1 - 4 

The  error  criterion  was  varied  from  10  to  10  to 
observe  the  effects  on  the  solution  and  computer  requirements. 
For  Crank-Nicolson  the  steady  state  solution  value  remained 
the  same  regardless  of  the  error  criterion.  This  demonstrated 
that  the  extra  time  required  to  satisfy  a tighter  error  cri- 
terion at  each  iteration  was  not  necessary  to  reach  a satis- 
factory steady  state  solution.  As  was  expected,  the  computer 
processing  time  did  increase  significantly  as  the  error 

criterion  was  made  smaller.  An  unusual  result  occurred  when 

-4 

the  error  criterion  was  set  to  10  . This  particular  value 

resulted  in  a processing  time  less  than  that  required  for 

10  1.  It  would  appear  that  this  might  have  been  the  result 

of  two  contributing  factors.  With  the  closeness  of  the  error, 

each  time  an  iteration  was  performed,  it  was  using  a better 

solution,  and  each  new  time  step  also  had  a better  solution 

value.  Although  there  were  many  more  time  steps  required 

for  this  error  criterion,  the  steady  state  value  was  rapidly 

approached  because  fewer  iterations  were  required  at  each  new 

time.  These  results  are  shown  in  Table  II.  These  results 

- 1 - 4 

for  an  error  criterion  between  10  and  10  do  not  imply 
that  this  would  be  the  result  obtained  for  all  problems  of 
this  type. 
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D.  EFFECT  OF  INCREASE  IN  TIME  STEP  SIZE 


In  this  method,  the  time  step  can  be  increased  or  not 
depending  on  the  solution.  If  a solution  is  not  obtained 
with  a particular  time  step,  that  same  time  step  is  reduced; 
and  the  routine  attempts  to  attain  solution  with  the  smaller 
time  step.  This  continues  until  either  a satisfactory  solu- 
tion value  is  obtained  or  the  built  in  default  value  is 
reached  which  terminates  the  routine.  When  a solution  is 
attained,  the  routine  compares  the  number  of  iterations 
required  to  reach  that  value  to  that  of  the  previous  time 
step  and  to  a specified  number  of  iterations.  If  the  current 
number  of  iterations  is  less  than  either  of  those  numbers, 
the  program  increases  the  time  step  by  an  input  constant 
value . 

The  change  made  to  the  initial  time  step  as  the  steady 
state  was  approached  was  the  critical  value  if  the  solution 
was  to  converge  to  a steady  state.  It  appeared  that  Crank- 
Nicolson  was  particularly  sensitive  to  the  amount  the  time 
step  was  increased  as  the  solution  was  approached  (Figure  9) . 
The  method  was  unable  to  recover  if,  when  nearing  the  knee  of 
the  curve,  the  size  of  the  increase  was  such  that  the  steady 
state  value  was  exceeded.  Once  the  solution  was  passed,  the 
results  indicated  a diverging  oscillation  about  the  steady 
state  value. 


Several  values  of  the  initial  time  step  were  utilized  in 
all  grids  to  attempt  to  locate  the  steady  state.  The  results 
of  these  trials  are  provided  in  Table  III.  It  is  noteworthy 


that  if  the  increase  was  greater  than  1.2  times  the  initial 
step  size,  the  steady  state  would  never  be  achieved.  Whether 
this  would  be  the  case  for  all  grids  is  not  obvious.  This  is 
particularly  evident  in  view  of  the  results  obtained  by 
Olsen  [6],  who  utilized  an  increase  of  1.5  times  the  initial 
step  for  a thirty-eight  node  grid.  It  would,  therefore,  seem 
that  a trial  and  error  approach  would  be  necessary  until  an 
increase  in  step  size  resulted  in  a steady  state  solution. 


IV.  DVOGER  (GEAR)  METHOD  OF  SOLUTION 

A.  DESCRIPTION 

The  DVOGER  (Gear)  routine  is  an  IMSL  library  routine 
which  integrates  a system  of  explicit  first  order  differ- 
ential equations.  Within  this  library  routine  is  the  capa- 
bility of  solving  both  stiff  and  nonstiff  systems  by  selection 
of  an  indicator.  In  solving  the  nonstiff  system,  the  Adams 
predictor-corrector  is  utilized.  For  the  stiff  system.  Gear’s 
predictor-corrector  is  used.  Gear’s  method  computes  the 
Jacobian  for  the  system  of  ordinary  differential  equations 
in  order  to  optimize  the  time  step  for  each  integration. 

B.  EFFECT  OF  DEGREES  OF  FREEDOM 

Gear's  method  placed  a much  greater  demand  on  the  computer 
in  terms  of  processing  time  and  storage  than  the  other  methods 
considered.  This  was  primarily  due  to  three  causes:  1)  the 

calculation  of  the  Jacobian,  2)  the  transformation  of  Equation 
(3)  to  explicit  form,  and  3)  the  initial  absolute  value  of 
the  solution.  The  first  two  items  require  a substantial 
increase  in  both  CPU  time  and  computer  storage.  Several 
variances  were  implemented  in  this  program  to  determine  if 
better  use  could  be  made  of  the  method  to  make  it  competitive 
with  Crank-Nicolson  or  the  Implicit  Gear  methods. 

When  the  system  was  considered  nonstiff,  the  storage 
requirements  decreased  since  the  Jacobian  was  not  calculated; 
but,  the  processing  time  increased  by  one  order  of  magnitude 
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when  compared  to  the  stiff  system  treatment  due  to  the  small 
time  steps  taken  to  solution.  The  routine  calls  for  an 
initial  absolute  solution  value  of  one  to  be  supplied.  This, 
however,  adversely  affects  the  progress  of  the  problem  since 
this  value  of  one  is  initially  compared  to  values  of  lO1^. 

This  limits  the  time  step  taken  by  the  routine  initially  and 
continues  to  affect  the  progress  until  the  steady  state  solu- 
tion is  neared.  When  an  initial  value  of  10*4  was  utilized, 
the  processing  time  was  decreased  by  forty  percent.  This  was 
due  to  the  fact  that  larger  time  steps  were  allowed  from  the 
outset,  since  the  previous  solution  value  (in  this  case  the 
initial  value)  was  comparable  to  the  value  obtained  in  the 
first  calculation. 

At  best,  when  utilizing  Gear's  method,  the  processing 
times  were  two  orders  of  magnitude  greater  than  that  required 
for  the  Implicit  Gear  method  and  twenty  times  the  Crank- 
Nicolson  processing  times  for  the  132  DOF  system,  as  shown  in 
Table  I.  Core  requirements  were  also  increased  significantly, 
particularly  at  one  hundred  thirty-two  degrees  of  freedom. 

For  this  grid,  DVOGER  (Gear)  was  approximately  three  times 
Crank-Nicolson  and  the  Implicit  Gear  core  requirements. 

Gear's  method  provided  the  same  steady  state,  solution 
values  as  the  other  two  methods,  and  the  transient  curves 
were  very  similar  to  those  obtained  in  the  Implicit  Gear 
method  (Figures  10-15). 
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C.  EFFECT  OF  ERROR  CRITERION 

In  Gear's  method,  the  time  step  size  is  adjusted  so  that 
the  single  step  error  estimate  divided  by  the  previous  maximum 
solution  value  is  less  than  the  error  criterion  in  the 
Euclidean  norm.  The  single  step  error  estimate  is  a multiple 
of  the  difference  between  the  predicted  and  corrected  values 
of  the  variable. 

The  error  criterion.  Table  II,  was  varied  for  this  method 
as  it  was  for  the  others.  The  results  were  the  same;  as  the 
criterion  was  tightened,  CPU  time  increased.  In  this  case, 
however,  the  time  became  excessive  very  rapidly,  more  than  an 


V.  IMPLICIT  GEAR  METHOD  OF  SOLUTION 


A.  DESCRIPTION 

This  numerical  method  is  particularly  useful  for  the 
solution  of  large,  sparse  systems  of  implicit,  stiff  differ- 
ential equations  [4].  In  contrast  to  Gear's  method,  the 
Implicit  Gear  method  eliminates  the  need  of  an  exact  NxN 
Jacobian,  but  rather  requires  only  an  exact  Nx7  Jacobian. 

This  is  due  to  the  fact  that  Equations  (3)  are  handled  directly, 
thereby  retaining  the  sparseness  of  the  matrix. 

Gear's  predictor-corrector  method  and  the  compact  matrix 
makes  use  of  storage  because  the  implicit  Equations  (3)  are 
handled  directly.  This  results  in  very  efficient  use  of  the 
computer  both  in  terms  of  storage  and  processing  time  require- 
ments . 

The  user  is  required  to  provide  a subroutine  that  evalu- 
ates the  system  of  equations  being  investigated,  as  well  as, 
for  efficiency,  a subroutine  to  evaluate  the  Jacobian  of  the 
system  of  equations. 

This  program  is  a modification  by  Franke  [4]  to  the 
routine  DFASUB  [7].  One  of  the  major  changes  to  the  routine 
is  in  the  treatment  of  the  erTor.  In  the  Implicit  Gear 
method,  instead  of  using  the  Euclidean  norm,  the  root  mean 
square  norm  is  utilized.  This  is  no  more  than  the  Euclidean 
norm  divided  by  the  square  root  of  the  number  of  components. 

In  addition,  the  maximum  value  of  the  component  is  updated 
before  the  norm  of  the  relative  error  is  computed. 
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B.  EFFECT  OF  DEGREES  OF  FREEDOM 


The  steady  state  values  obtained  in  this  method  were  the 
same  as  those  obtained  in  Crank-Nicolson  and  DVOGER  (Gear) 
as  shown  in  Figures  11-18.  This  result  does  not  imply  that 
these  methods  will  always  give  identical  steady  state  solu- 
tions. Certainly  they  do  give  different  transient  solutions. 
The  major  advantages  of  this  method  are  1)  the  computer 
processing  time  was  reduced  as  the  grid  size  became  finer 
and  2)  the  storage  requirements  are  slightly  greater  than 
that  required  in  the  Crank-Nicolson  method  due  to  the  compu- 
tation of  the  Jacobian  and  the  storage  of  up  to  seven  past 
solutions  which  control  the  size  of  the  time  step.  Implicit 
Gear  requires  about  Nx2S  more  storage  locations  than  Crank- 
Nicolson.  This  amounts  to  about  only  thirteen  thousand  bytes 
for  the  one  hundred  thirty- two  degrees  of  freedom  system 
(single  precision).  The  magnitudes  of  the  increase  in  pro- 
cessing time  dropped  dramatically  in  this  method.  For  the 
same  initial  disturbance,  the  time  increased  eighty  percent 
for  the  one  hundred  thirty- two  node  grid  and  twenty- four 
percent  for  the  seventy-two  degrees  of  freedom.  In  addition, 
Implicit  Gear  gave  a more  accurate  transient  solution.  The 
results  of  this  method  are  shown  in  Table  I and  comparison 
to  the  other  methods  will  be  made  in  Chapter  VI. 

C.  EFFECT  OF  ERROR  CRITERION 

In  Implicit  Gear,  the  error  criterion  is  defined  as  in 
DVOGER  (Gear)  except  that  the  root  mean  square  of  the  Euclidean 
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norm  and  the  updated  maximum  solution  value  are  used  in  the 
computation. 

The  error  criterion  was  varied  substantially  to  determine 

what  effect  it  had  on  computer  storage  and  time  requirements, 

as  well  as  its  effect  on  the  accuracy  of  solution.  As  shown 

in  Table  II  and  Figures  19-20,  a wide  range  of  convergence, 

-4 

1.0  to  10  , was  investigated  with  interesting  results. 

Expectedly,  the  processing  time  did  increase  as  a tighter 

error  criterion  was  required.  However,  the  criterion  had 

little  effect  on  the  track  through  the  transient  solution  for 

values  of  less  than  10*.  When  using  an  error  of  one-half 

and  one,  the  transient  solution  varied  substantially.  The 

same  steady  state  solution  was  obtained,  although  at  a later 

point  in  problem  time  as  shown  in  Figure  19.  For  the  stiffer 

error  criterion  requirements,  the  processing  time  necessary 

- 4 

is  doubled  by  utilizing  a criterion  of  10  instead  of  a value 
greater  than  10"*.  This,  certainly,  did  not  appear  to  warrant 
the  closer  tolerance  level. 
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VI . RESULTS 


A.  COMPARISON  OF  TIME  AND  STORAGE  REQUIREMENTS 

With  all  three  methods  yielding  the  same  steady  state 
solution  (Figures  10-18),  the  comparison  of  methods  became 
one  of  time  and  storage  requirements  placed  on  the  computer 
rather  than  one  based  on  which  provided  the  best  solution. 

In  the  transient  stage,  Gear  tracked  slightly  better  than 
Implicit  Gear,  but  both  Gear  and  Implicit  Gear  tracked  better 
than  Crank-Nicolson  (See  Figures  10-18).  For  the  three  systems 
of  equations  utilized,  the  Implicit  Gear  method  performed 
significantly  better  than  either  Crank-Nicolson  or  DVOGER 
(Gear)  in  CPU  time  and  was  comparable  to  Crank-Nicolson  and 
superior  to  Gear  in  storage  requirements  (Table  I) . This 
becomes  more  significant  as  the  number  of  degrees  of  freedom 
increase.  The  Implicit  Gear  is  less  sensitive  to  the  increase 
in  size  of  the  system  of  equations  as  shown  in  Figure  8. 

Comparing  the  forty-five  and  the  one  hundred  thirty-two 
degrees  of  freedom,  the  core  requirements  increased  slightly 
and  the  time  doubled  for  the  Implicit  Gear  method.  In  con- 
trast, for  Crank-Nicolson,  the  time  increased  by  six  times 
and  the  core  by  100K  bytes;  and,  for  DVOGER  (Gear),  there 
was  a ten  times  and  500K  byte  increase  in  time  and  storage 
respectively. 
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B.  DEGREE  OF  FREEDOM  EFFECT  ON  SOLUTION 

Varying  the  degrees  of  freedom  resulted  in  a better 
solution  value  as  shown  in  Figures  5-7  and  21-26.  There  was 
a four  and  one-half  percent  better  solution  obtained  for  the 
one  hundred  thirty-two  degrees  of  freedom  compared  to  the 
forty-five  degrees  of  freedom.  This  better  solution  is  very 
costly  both  in  terms  of  computer  core  requirements  and  CPU 
time  (Figure  8) . 

C.  ERROR  CRITERION 

Error  criterion  was  varied  for  all  methods  to  determine 
the  effect  on  solution  and  computer.  As  might  be  expected, 
the  computer  processing  time  increased  for  all  methods  as 
more  rigid  tolerance  was  imposed  (Table  II).  However,  although 
the  criterion  was  made  very  close,  there  was  no  appreciable 
effect  on  the  transient  solution  until  the  error  criterion 
was  greater  than  10  1 and  no  effect  on  the  steady  state 
solution.  This  is  significant  in  that,  as  stated  above,  the 
processing  time  increases  greatly  as  the  criterion  is  tight- 
ened. This  shows  that  for  the  particular  problem  considered 
here,  some  close  error  criteria  yield  no  advantage,  only  the 
disadvantage  of  requiring  more  time  in  the  computer. 

D.  INITIAL  DISTURBANCES 

The  input  of  various  initial  disturbances  (central,  skewed 
at  the  core- reflector  interface  and  uniform  throughout  the 
core)  yielded  the  same  steady  state  solution  value.  The 
transient  solution  varied  substantially  due  to  the  input 
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position  of  the  disturbance.  In  most  cases,  due  to  the 
magnitude  and  scope  of  the  initial  disturbance,  the  uniform 
disturbance  required  the  least  processing  time  to  reach  the 
steady  state  value. 

E.  COMPARISON  OF  SOLUTION  TRACKING 

Two  nodal  points  were  investigated  graphically  for  each 
disturbance,  integration  method  and  DOF,  as  shown  in  Figures 
5-7,  10-18  and  21-26.  At  these  test  points,  the  methods 
performed  similarly  in  all  cases  in  both  the  transient  and 
steady  state  solutions. 
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VII.  TEST  PROBLEMS 


A.  PHYSICAL  MODEL 

A cylindrical  reactor  with  the  dimensions  and  properties 
given  in  Table  IV  was  used  as  a model.  A radial  slice  of 
this  model  was  discretized  by  various  finite  element  grids 
as  shown  in  Figures  1-3. 

B.  COMPUTER  PROCESSING  CONSIDERATIONS 

The  programs  presented  in  this  thesis.  Appendix  A,  were 
written  in  the  FORTRAN  IV  language  and  all  computer  runs. 

Table  V,  processed  on  the  IBM  360/67  computer  using  the 
FORTRAN  ' G * compiler.  It  is  expected  that  the  FORTRAN  'H' 
compiler  would  have  supplied  results  similar  to  those  obtained 
with  the  FORTRAN  'G'  compiler.  This  was  not  done  because  the 
FORTRAN  'H'  compiler  requires  350K  bytes  as  a minimum  core 
requirement.  The  majority  of  the  runs  performed  in  this 
thesis  required  significantly  less  than  300K  storage.  Single 
precision  (six  to  seven  significant  digits)  was  used  through- 
out this  research.  As  has  been  shown  [6] , double  precision 
solutions  are  at  variance  by  less  than  0.01  percent  with 
single  precision  solutions. 

C.  PROBLEM  ANALYSIS 

The  techniques  used  by  Olsen  [6] , modified  as  described, 
and  the  Implicit  Gear  method  [4]  were  utilized  to  solve  the 
problem  delineated  in  Ref.  [5] . The  problem  was  approached 
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using  finite  element  discretizations  with  forty-five,  seventy- 
two  and  one  hundred  thirty- two  degrees  of  freedom  and  providing 
an  initial  disturbance.  Three  initial  disturbances  were  pro- 
vided as  follows:  central  at  the  core  center  (R  = 0 cm., 

Z * 0 cm.);  uniform  throughout  the  core;  and,  a skewed  dis- 
turbance at  the  core-reflector  interface  (R  = 60  cm.,  Z = 0 
cm. ) . 

D.  PROGRAM  USAGE 

In  order  to  properly  utilize  the  routines  presented  in 
this  thesis,  several  points  must  be  considered. 


In  the  Crank-Nicolson  routine,  there  was  difficulty  in 
obtaining  a steady  state  solution.  The  ultimate  cause  was 
the  size  with  which  the  previous  time  step  increases.  This 
requires  a trial  and  error  approach  for  the  particular  grid 
size.  For  this  project,  an  increase  of  ten  or  twenty  percent 
yields  satisfactory  results;  but,  anything  larger  results  in 
a divergent  oscillation  about  the  steady  state  value. 

In  all  three  methods,  the  dimensioning  should  be  set  by 
the  particular  problem  being  solved,  i.e.,  determined  by  the 
degrees  of  freedom.  If  this  is  done,  most  efficient  use  will 
be  made  of  the  computer,  and  the  user  will  experience  better 
turn  around  time  on  the  equipment. 

In  the  use  of  the  Implicit  Gear  method,  the  data  cards 
for  nodal  neighbor  connectivity  must  have  the  contributing 
nodes  listed  consecutively , with  the  zeros,  used  for  nodes 
not  having  six  neighbors,  being  last  on  the  cards.  For  example, 
in  the  one  hundred  thirty- two  degrees  of  freedom  system 
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(Figure  3) , the  nodal  neighbor  connectivity  for  node  122 

would  be  written  as 

122  123  11  0 0 00. 

The  data  deck  arrangement  for  all  three  methods  is  given  in 
Table  VI . 
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VIII.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

It  has  been  shown  in  this  research  project  that  of  the 
three  methods  investigated  in  the  solution  of  the  test  problem 
that  the  Implicit  Gear  performed  in  a superior  manner  in  CPU 
time  requirements  and  was  comparable  to  the  storage  require- 
ments of  Crank-Nicolson.  In  addition,  Implicit  Gear  was  the 
least  sensitive  to  the  change  in  the  degrees  of  freedom. 

The  method  used  had  no  effect  on  the  steady  state  solution 
value  obtained  as  all  three  yielded  the  same  results  fo.r  the 
same  DOF.  In  the  transient  solution,  Implicit  Gear  conformed 
very  closely  to  DVOGER  (Gear).  As  shown  in  Figures  10-15,  as 
the  DOF  increased,  the  transient  solutions  moved  in  the  direc- 
tion of  the  transient  solution  of  Gear's  method. 

As  the  DOF  in  the  discretized  finite  element  was  increased, 
a better  solution  was  obtained.  This,  however,  was  counter- 
acted by  the  fact  that  for  the  small  increase  in  accuracy 
obtained  by  the  finer  mesh,  the  time  and  storage  need  of  the 
computer  increased  significantly. 

A very  important  result  of  this  project  was  the  effect  of 

varying  the  error  requirements  of  the  problem.  This  yielded 

essentially  no  change  in  the  accuracy  of  the  transient  solu- 

- 1 - 4 

tion  for  a range  of  10  to  10  and  no  variance  in  the  steady 

- 4 

state  solution  for  a range  from  1.0  to  10  . In  using  these 

methods  of  integration,  the  error  criterion  selected  would 
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be  a function  of  the  solution  desired.  If  the  user  is  only 
concerned  with  the  steady  state  solution,  the  tolerance  could 
be  relaxed  to  1.0.  However,  if  the  transient  solution  is 
needed,  an  error  of  10  * appears  to  be  the  least  rigid  value 
which  will  still  provide  a satisfactory  track  to  steady  state. 


B.  RECOMMENDATIONS 

It  would  be  of  value  to  investigate  further  the  Implicit 
Gear  method  through  the  use  of  other  problems  solving  a non- 
linear system  of  ordinary  differential  equations.  In  addition, 
larger  degrees  of  freedom  should  be  attempted  to  evaluate  any 
limitations  that  may  exist  in  this  method  of  integration. 


TABLE  I 


COMPARISON  OF  SOLUTION  METHODS 


METHOD 

CORE 

(K) 

PROCESSING  TIME* 
Central  Skewed 

(min. ) 
Uniform 

45 

Nodes 

CRANKO 

158 

1.67 

1.73 

3.57 

DVOGER  (GEAR) 

148 

17.8 

19.5 

12.5 

IMPLICIT  GEAR 

124 

1.05 

1.00 

1.20 

72 

Nodes 

CRANKO 

184 

3.07 

5.49 

2.38 

DVOGER  (GEAR) 

244 

89.5 

114.1 

59.5 

IMPLICIT  GEAR 

124 

1.3 

1.30 

1.25 

132 

Nodes 

CRANKO 

252 

9.1 

9.1 

9.1 

DVOGER  (GEAR) 

610 

>400 

>400 

>400 

IMPLICIT  GEAR 

124 

1.9 

2.0 

2.0 

I I 

* Processing  time  is  that  required  to  reach  a steady  state 
solution. 
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TABLE  II 


1 


EFFECTS 

OF  SOLUTION  ERROR 

CRITERION 

SOLUTION  ERROR 
CRITERION 

45  Nodes 

PROCESSING  TIME 
(rain. ) 

TIME  AT  SOLUTION 
(sec. ) 

CRANK -NI COLSON 

0.1 

1.8 

3.69  x 10‘5 

0.01 

5.5 

3.70  x 10'5 

0.001 

14.2 

3.57  x 10'5 

0.0001 

0.83 

3.87  x 10*5 

IMPLICIT  GEAR 

1.0 

0.99 

7.56  x 10'4 

0.5 

0.92 

1.48  x 10'2 

0.1 

1.05 

7.06  x 10'5 

0.01 

1.20 

7.64  x 10'5 

0.  001 

1.45 

6.36  x 10'5 

0.0001 

1.58 

6.44  x 10'5 

DVOGER  (GEAR) 

0.1 

18.6 

5.64  x 10"5 

0.01 

29.5 

5.79  x 10~5 

0.001 

78.8 

5.55  x 10'5 

0.0001 

>300 

_ 
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TABLE  III 


EFFECTS  OF 

VARYING  THE  INCREASE 
IN  CRANK- N I COLSON 

OF  THE 
METHOD 

INITIAL  TIME  STEP 

INCREASE 

CORE 

PROCESSING 
(rain. ) 

TIME 

SOLUTION  TIME 
(sec. ) 

45  NODES 

1.1 

158K 

4.99 

3.18  x 10‘5 

1.2 

158K 

1.65 

2.69  x 10"5 

1.5 

184K 

7.42 

OSCILLATING 

72  NODES 

1.1 

252K 

3.5 

4.38  x 10~5 

1.2 

184K 

3.05 

3.23  x 10'5 

1.3 

252K 

3.5 

OSCILLATING 

1.4 

252K 

9.1 

OSCILLATING 

1.5 

184K 

6.1 

OSCILLATING 

132  NODES 

1.1 

252K 

9.1 

5.58  x 10'5 

1.2 

252K 

8.4 

8.03  x 10‘5 

1.5 

184K 

2.6 

OSCILLATING 

TABLE  IV 


PHYSICAL  CONSTANTS 


SYMBOL 


COMPUTER 

SYMBOL 


DEFINITION 


VALUE 


C-N,  D / IG 


R 

R 

Total  radius 

90  cm. 

Rc 

R 

Core  radius 

60  cm. 

Hc 

Z 

Core  height 

160  cm. 

H 

z 

Total  height 

220  cm. 

V 

V/VELOCT 

Neutron  velocity 

7 

4.8x10  cm/sec 

Dc 

D/DSUBF 

Neutron  diffusion 
coefficient  (core) 

0.913  cm. 

Dr 

D/DSUBC 

Neutron  diffusion 

coefficient 

(reflector) 

1.20  cm. 

Zac 

SGA/SGMAAF 

Neutron  absorption 

cross-section 

(core) 

0.01401  cm-1 

Zar 

SGA/SGMAAC 

Neutron  absorption 

cross-section 

(reflector) 

0.008  cm"1 

V 

ZNU 

Number  of  neutrons 
per  fission 

2.54 

Zfc 

SGF/SGMAF 

Neutron  fission 

cross-section 

(core) 

0.008  cm'1 

Efr 

SGF/  - 

Neutron  fission 
cross-section 

0.0  cm  1 

e 

FISFAC 

Fission  energy 

7.652  x 10"12  cal/fis 

Ka/v 

HBAR 

Modified  convection 
heat  transfer 
coefficient 

0.0632  cal/cm3  -sec-°C 

a 

ALPHA 

Reactivity  temper- 
ature coefficient 

1 x 10"5  / °C 
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TABLE  V 


1 


RUN 

METHOD 

NODES 

DISTURBANCE 

COVERGENCE 

TIME  STEP 
CHANGE 

1 

CRANK- NI COLSON 

45 

CENTRAL 

0.1 

1.1 

2 

1.2 

3 

1.5 

4 

0.01 

1.2 

5 

0.001 

6 

0.0001 

7 

UNIFORM 

0.1 

8 

SKEWED 

9 

72 

CENTRAL 

0.1 

1.0 

10 

1.1 

11 

1.2 

12 

1.3 

13 

1.4 

14 

1.5 

15 

UNIFORM 

1.2 

16 

SKEWED 

17 

132 

CENTRAL 

0.1 

1.1 

18 

1.2 

19 

1.5 

20 

UNIFORM 

1.2 

21 

SKEWED 

22 

DVOGER 

MTH=0 

45 

CENTRAL 

0.1 

N.  A. 

23 

DVOGER 

MTH-1 

YMAX-1 

0.1 

45 

TABLE  V (Continued) 


RUN 

METHOD 

NODES 

DISTURBANCE 

CONVERGENCE 

24 

0.01 

25 

0.001 

26 

DVOGER 

MTH*1 

YMAX*1 

45 

CENTRAL 

0.0001 

27 

UNIFORM 

0.1 

28 

SKEWED 

29 

72 

CENTRAL 

30 

UNIFORM 

31 

SKEWED 

32 

132 

CENTRAL 

33 

UNIFORM 

34 

SKEWED 

35 

IMPLICIT 

GEAR  45 

CENTRAL 

0.1 

36 

0.01 

37 

0.001 

38 

0.0001 

39 

UNIFORM 

0.1 

40 

SKEWED 

41 

72 

CENTRAL 

0.1 

42 

UNIFORM 

43 

SKEWED 

44 

132 

CENTRAL 

0.1 

45 

UNIFORM 

46 

SKEWED 

47 

45 

CENTRAL 

0.5 

TIME  STEP 
CHANGE 
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TABLE  V (Continued) 


RUN 

METHOD 

NODES 

DISTURBANCE 

CONVERGENCE 

48 

1.0 

49 

DVOGER 

mth-i  14 

YMAX»10ah 

45 

CENTRAL 

0.1 

TIME  STEP 
CHANGE 
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TABLE  VI 


DATA  DECK  ARRANGEMENT 

CRANK- NICOLSON 
Title 

NUMEL ,NUPBP , NUMSNP , NFULEL 
List  of  outer  boundary  points 
MTH,MAXDER,NCOUNT 

ZNU.FISFAC, HBAR , EPSVAL , ERRVAL , AFUEL 
TO ,H ,TF ,HMIN ,HMAX 

V - neutron  velocity 
D - diffusion  length 

SGA  - absorption  cross-section 

?GF  - fission  cross-section 

AIIHA  - reactivity  temperature  coefficient 

PSIIV  - initial  disturbance  flux 

System  nodal  points  with  axial  and  radial  position 
Element  nodal  connectivity 
Nodal  neighbor  connectivity 

DVOGER 

Title 

NUMEL , NUPBP , NUMSNP , NFULEL 
List  of  outer  boundary  points 
MTH,MAXDER,NCOUNT 

ZNU , F I SF AC , HBAR , EPSVAL , ERRVAL , AFUEL 
TO,H,TF,HMIN,HMAX 

V 
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TABLE  VI  (Continued) 

D 

SGA 

SGF 

ALPHA 

PSIIV 

System  nodal  points  with  radial  and  axial  position 
Element  nodal  connectivity 


IMPLICIT  GEAR 


NCORE 

List  of  core  elements 


with  uniform  initial  disturbance 


Title 


NUMEL , NBP , NUMSNP , NFULEL , NELROW ,MXE VAL , NCOUNT , NCMPK 
NDE.NL 

List  of  outer  boundary  points 

VELOCT , DSUBF , DSUBC , SGMAAF , SGMAAC , SGMAF , FISFAC ,HBAR 
ZNU .ALPHA , RMSEPS , AFUEL .TSTART .TEND , PINITV .SCALE 
HMIN.HMAX 

MTH , MAXDER , NPROB , NTYPE , JPLOT , NPLOT , I RUN 
Nodal  neighbor  connectivity 

System  nodal  points  with  radial  and  axial  position 
Element  nodal  connectivity 
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Z cm. 


Figure  1.  Reactor  Model  with  45  Nodes 
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Figure  3.  Reactor  Model  with  132  Nodes 
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60  cm 


Time  Dependent  Neutron  Flux  for  All  DOF  Using  Crank-Nicolson 
Method  for  a Central  Disturbance  (R  = 0 cm.,  Z = 0 cm.) 


Time  Dependent  Neutron  Flux  for  All  DOF  Using  Crank-Nicolson  Method 
for  a Skewed  Disturbance  (R  = 60  cm.,  Z = 0 cm.) 


Effect  of  Time  Step  Change  in  Crank-Nicolson  on  Achieving  Steady  State. 


Figure  10.  Time  Dependent  Neutron  Flux  for  All  Methods  with  45  DOF 
for  a Central  Disturbance  (R  = 0 cm.,  Z = 0 cm.). 


Dependent  Neutron  Flux  for  All  Methods  with  45  DOF 
Skewed  Disturbance  (R  = 60  cm..  Z * 0 cm.l 


Dependent  Neutron  Flux  for  All  Methods  with  45  DOF 
Uniform  Disturbance  Throughout  Core. 


ZOC 


Figure  13.  Time  Dependent  Neutron  Flux  for  All  Methods 
for  a Central  Disturbance  (R  = 0 cm..  Z = n 


Z0C.<W 


Dependent  Neutron  Flux  for  All  Methods  with  72  DOF 
Skewed  Disturbance  (R  = 60  cm.  , Z = 0 cm.'). 


Z6e.*i 


Dependent  Neutron  Flux  for  All  Methods  with  72  DOF 
Uniform  Disturbance  Throughout  Core. 


ICQ**- 


( pO  OA1 


Figure  19.  Comparison  of  Error  Criterion  for  Implicit  Gear  Method. 


23.  Time  Dependent  Neutron  Flux  for  All  DOF  Using  Implicit  Gear 
Method  for  a Uniform  Disturbance  Throughout  the  Core. 


jure  24.  Time  Dependent  Neutron  Flux  for  45  and  72  DOF  Using  the 
DVOGER  (Gear)  Method  for  a Central  Disturbance 


Figure  25.  Time  Dependent  Neutron  Flux  for  45  and  72  DOF  Using  the 
DVOGER  (Gear)  Method  for  a Skewed  Disturbance 


Time  Dependent  Neutron  Flux  for  45  and  72  DOF  Using  the 
DVOGER  (Gear)  Method  for  a Uniform  Disturbance  Throughout 
the  Core. 
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SUBROUTINE  PEPSI  ( NUMEL , NUMSNP, R I ,R2 , R3 , AR EA, ELNOO 
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FINITE  ELEMENT  SOLUTION  OF  NONLINEAR  REACTOR  DYNAMICS 
IN  TWO  DIMENSIONAL  SPACE  (COMPACT  VERSION! 
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oooo  o on  noon 


—DATA  GEMERATOR  TO  PROVIDE  SYSTEM  MODAL  POINT 
RADIAL  AMO  AXIAL  POSITION 


20 

10 


IMPLICIT  REAL  * 4 ( A-H.O-Z) , INTEGER  ♦ 4 II-N) 
INTEGER  * 4 Sf SNODf ELNOD 

DIMENSION  RC 1 53  I , Z(15  0) , X( 150 ) , Y ( 150) , SYSN0D( 150) , 


1 

100 

200 

300 

400 


500 

501 


60 

700 


600 

601 


1 ELN  JD ( 220 « 3 ) 

READ  IN  GRID  VALUES  AND  NJMBER  OF  NODES 

READ! 5* 100)  NH»  MV 
FORMAT (215) 

WRI TE  ( 6 * 200) 

FORMAT  Cl*  ) 

WRITE ( 6 * 300)  NH.NV 

FORMAT ( 2X*  *NH  *'tI5,5X,'NV  =',I5) 

READ!  5»  40  0)  ( X(  I)  ♦ 1*1  , NH ) 

READ(  5*  400)  ( Y( J ) , J*1 * NV) 

FORMAT ( 15F5*1 ) 

WRITE (6  »500)  (X( I ) t I = 1»NH) 

WRITE (6*501)  (Y( J)  ,J  = i*NV) 

FORMAT ( //2X» • R ' *5X*14F8«3) 

FORMAT (//2X»'Z,*5X*14F8.3) 

NVH=NV*NH 

NN*NV-1 

K=0 

DO  10  1*1, NH 
DO  20  J*1,NN 
K=K+ 1 

SYSN30( K )=K 
R ( K) *X ( I) 

Z(K)*Y( J) 

CONTINUE 
CONTINUE 
K*NVH-NH 
DO  60  I *1 » NH 
K*K+1 

SYSNOD ( K ) =K 
R(K)=X( I) 

Z(K)=Y( NV) 

CONTINUE 
WRITE (6,700) 

FORMAT( /////5  < * * NODE • ,11X , • R« , 1 IX , • Z* ) 
WRITE (6, 600)  (I,R( I), Z( I) ,I*1,NVH) 

WRITE (7, 601)  (I,R(I)fZ(I) ,I=1,NVH) 

FORMA T( /5X,I5,2F15.5) 

FORMAT (5X*I5*2F15*5) 


40 

30 


COMPUTE  NODAL  POINT  CONNECTIVITY  OF  ELEMENTS 

LL*  1 

K*1 

MM*0 

NHH*NH- 1 

NVV*NV-2 

DO  30  1*1, NHH 

DO  40  J*1*NVV 

ELNOD (K  * 1 J =LL 

E LNOD ( K+l • 1 ) = wL 

ELNOO(  K *2)  *MM*'NV 

ELNOD (K+l • 2 ) = LL+NV 

ELNOD (K,3)*ELN0D(K+1,2) 

ELN0D(K+l,3)*LL*l 

LL*LL«-i 

K*K+2 

MM*MM*1 

CONTINUE 

LL*LL*1 

MM=MM+1 

CONTINUE 
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^ 


70 

750 

800 

950 


50 

900 

901 


DO  70  I * 1 • NHH 
ELNODtK  tl i *NN 
ELN0D(K+1.1)=NM 
ELN00(K,2)»NN+NV-l 
ELNQD(K+1 »2 )SLL+NV 
ELN0D1 K »3J=ELN0D ( K+l »2 ) 
ELN0D(K+it3»**tm-NV 
K =K+2 

NN=*NN+NV- 1 

LL*LL+1 

MM=MM*1 

NEL=K-1 

CONTINUE 

WRITE ( 5 *750) 

FORMAT! • 1' ) 

WRITE(6,800)NEL 

FORMA T(5X»  * NUM3ER  OF  ELEMENTS  **,I4) 

WR ITE ( 6 *950 ) 

FORMAT(/////10X, 'CONNECTIVITY  MATRIX* ) 


DO  50  I * i • NEL 
WR ITE ( 6 1 900 ) 1 1 ( ELNOD (It 
WRITE ( 7 *901)  I » ( ELNOD ( I 
CONTINUE 
FORMAT !/4I 10 ) 
F0RMATC4I10) 
IF(NV.EQ.12)ST3P 
GO  TD  1 
END 


J=l*3) 
t J=l » 3) 
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non  o ooo 


NUCLEAR  PROPERTY  DATA  GENERATOR 


100 

200 

*00 

300 


10 


20 

30 


500 

600 


INTEGER  * 4 NEL . ELEM, ITYPE 

DIMENSION  V (220  I , ALPHA! 220 ) » SGF ( 220) ,S3A!220) ,0(220), 
1 ITYPEI220) ,ELEM(220) 

READ( 5, 100)  NEL 

IDENTIFY  CORE  AND  REFLECTOR  ELEMENTS 
READ(5. 200)  (ITYPE(I)  ,1=1, NEL) 

FORMAT (2110) 

FORMAT ( 1615) 

WRI TE ( 6*400) 

FORMAT! 1 1*, 2X,'ELEM», 6X, »D*  ,8X, • SGA • , 9X, • SGF* ,8X, 

1 'ALPHA' ,8X,' V' ) 

FORMAT! 2X, 14,5!  IPE12.4) ) 

ASSIGN  PROPERTY  VALUES  TO  ELEMENTS 

DO  30  1*1, NEL 

ELEM! I) *1 

V ( I )*4»8E+07 

ALPHAII )*1.0E-05 

IF! ITYPE(I) .EQ.O)GO  TO  10 

0! I )*1  .2 

SGA! I )=.008 

SGF ! 1 1*0 

GO  TO  20 

D ( I ) *0. 913 

SGA(I)=0.014 

SGF! I )=0.008 

WRITE! 6, 300)! 1,0! I) ,SGA!I ) ,SGF(I ), ALPHA! I) , V! I ) ) 
CONTINUE 

WRITE (7 ,500) ( V!  I), 1*1, NEL) 

WRITE!7,500) (D! I ) ,1*1, NEL) 

WRITE! 7,500) ($3A(  I ),  1*1, NEL) 

WRITE !7, 500) (SGF! I ), 1*1, NEL) 

WRITE! 7,600) ! AL PHA! 1 ) , I *1 ,NEL ) 

FORMAT ( 8G10*  5 ) 

FORMAT! 5 FI 5. 8) 

I F ( NEL. EQ* 64) GO  TO  1 
IF(NEL.EQ.112)G0  TO  1 
IF(NEL.EQ.220)ST0P 
END 


ooooo  o no  noon 


DATA  GENERATOR  TO  PROVIDE  NODAL  NEIGHBOR 
CONNECTIVITY 


1 

200 


INTEGER  * 4 NV,NVH,LCON,NN,MM, LL , JJ , MNOD , 1 1 ,KK, 
I I J » I K« IM,IN,NH,LI*  LJtLK 
DIMENSION  MN0D( 132  » 7) 

REAO  IN  INITIAL  DATA 
READ! 5* 200)  NV, N VH, LCON ,NH 


FORMAT (415  I 
COMPUTE 
M-NV+1 
N*NV-1 
NN*2*N 
MM*3*N 
1 1 *4*N 
KK»5*N 
IJ*6*N 
I K*7*N 
I M»8*N 
IN*9*N 
LK*10*N 
LI*NVH-NH 
LJ*LI+1 
LL*LI-NV+2 


MATRIX  COUNTERS 


10 

20 

30 

40 

50 

60 


CALCULATE 

MATRIX 


COMPACTED  NODAL  NEIGHBOR 


DO  10  1*1, 
MNOD( 1,1)* 
MNOD< 1,2)* 
MNOD( 1,31* 
MNOOI 1,4)* 
CONTINUE 
DO  20  1*2, 
MNOD( I, 5)* 
CONTINUE 
DO  30  I *NV 
MNOD(  1,6)* 
CONTINUE 
DO  40  I*M, 
MNOOI I, 7)* 
CONTINUE 
DO  50  1*1, 
MNOD (1,7)* 
CONTINUE 
DO  60  1*1, 
MNOD (I, 6)* 
CONTINUE 
MNOD( 1,5)* 
MNOD( N» 2) * 
MN0D(N,4)* 
MNOD( NV ,5 ) 
MNODINN ,2 ) 
MN0D(NN,4) 
MN0DINN+1, 
MNOOI NN+i • 
MNOD (MM, 2) 
MNOOI MM ,4) 
MNOD! MM+1 , 
MNOD( MM+1 , 
MNOOI LI +1 , 
MNOO(NVH,2 
MNOD! II ,2) 
MNOD! II ,4) 
MNOOI 1 1 +1, 
MNOD! 1 1 +1 , 
IFINVH.EQ. 
MNOD IKK, 2) 


NVH 

I 

1 + 1 
N+I 
NV+I 

NVH 

1-1 

, NVH 
I-N 

NVH 

I-NV 

NV 

0 

N 

0 

0 

LI  + l 
LI+2 
*0 

=LI  *-2 
*L  1*3 
5 ) *0 
7 )*0 
»LI+3 
*L  I *4 
5 ) *0 
7 ) * 3 
5)*0 
) =0 
*L  I 
*L  I *5 
5 )*D 
7 ) *3 
45)30 
=L  I *-5 


TO  15 


CONNECTIVITY 
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15 

70 

80 

85 


95 

300 


90 

100 

101 


MNOOl 

MNOO 

MNOO 

NNOD 

MNOO 

MNOO 


MN0D(KK,4)*Ll+6 
MNOO ( KK+i » 5 ) * D 
MNOOt KK+1  ,7)*D 
MNOOt IJ ,2)=LI+6 
MN0D(IJ,4)=LI+7 

IJ+1,5>0 

1 J +1 , 7 ) * 3 
IK,2)*LI+7 
IK,4)*LI+8 
IK+1,5)*0 

IK+1,7)*0 

IFtNVH.EQ. 72)30  TO  15 
MNOOt IM»2)=LI+8 
MNOOt  IM,4)*LI+9 
MNOOt IM+1,  5)*0 
MNOOt IM+1 • 7)  = 3 
MNOOt  LK, 2)  *LI +13 
MN0D(LK»4)=LI+li 
MNOOtLK+1 y 5) s0 
MNOOt  LK+1 *7 ) O 
MNOOt IN, 2>=LI+9 
MNOOt IN,4)*LI+13 
MNOOt IN+1,5)*0 
MNOOt IN+1 » 7 ) = 3 
00  70  I *LJ»  N\H 
MNOOt 1,6) >0 
MNOO  1 1 , 7 ) =0 
CONTINJE 
00  80  I *LL , Li 
MNOOt I, 3 )=0 
CONTINJE 
00  85  I *LL» LJ 
MNOOt I , 4) =6 
CONTINUE 
DO  95  I *L  J » NVM 
MNOOt 1,3) *N 
IF(I.EQ.NVH)G3  TO  95 
MNODt 1+1,4) =M')DD 11,3) 

N*N+NV-1 
CONTINUE 
MNOOt  LI  ,2)=NV-t 
WRITE (6,300) 

FORMA Tt  5X  »' HEXAGONAL  CONNECTIVITY',/ 
00  90  1*1, NVH 

WRITE  (6  »100M MNOO (I,J),J=l,LCON) 


WRITEt7,10i)tMN03(I,J) 
CONTINJE 

F ORMAT ( /9 ( 2X  » 1 3 ) ) 

FORMA Tt  9( 2X* 13) ) 
IFtNV.EQ.12)ST0P 
GO  TO  1 
END 


> J«l, LCON) 


• MODE ' 
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